#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxUtils.h"

void fittpc()
{
  //TPC pp
  //aliroot ../sep/outQ.root

  TH1D *hfit=(TH1D*)((TH1D*)gDirectory->Get("hhtpc_mpv"))->Clone("hfit");

  Double_t fPar[]={10, 4.4, 0.004, 2.26, 0.95};
  Int_t npar = 5;
  XGLUtils::ChisquareFit(hfit, AliTRDdEdxUtils::MeandEdxLogx, npar, fPar);
  
  TCanvas *cc=new TCanvas("cc");
  hfit->SetMarkerStyle(20);
  hfit->SetMarkerSize(1);
  hfit->Draw();

  TF1 *ff=new TF1("ff",AliTRDdEdxUtils::MeandEdxLogx, -0.5, 4.5, npar);
  ff->SetParameters(fPar);
  ff->SetLineColor(kRed);
  ff->Draw("same");

  cout<<endl;
  gSystem->Exec("date +%m%d%Y%H%M%S");
  for(Int_t ii=0; ii<npar; ii++){printf("par[%d]=%15f;\n", ii, ff->GetParameter(ii));} cout<<endl;cout<<endl;
  
  TFile *fout=new TFile("fittpc.root","recreate");
  hfit->Write();
  ff->Write();
  cc->Write();
  fout->Save();
  fout->Close();
}

void fittrd(const TString opt)
{

  //TRD pp/PbPb:

  const TGraphErrors *g1=(TGraphErrors*)gDirectory->Get("type1mpv");
  const TGraphErrors *g2=0x0;
  if(opt.Contains("pp")){
    g2 = (TGraphErrors*)gDirectory->Get("type2mpv");
  }
  else if(opt.Contains("PbPb")){
    g2 = (TGraphErrors*)gDirectory->Get("type4mpv");
  }
  else{
    printf("opt error!! %s\n", opt.Data());exit(1);
  }

  TGraphErrors * hcomb=new TGraphErrors;
  hcomb->SetName("hcomb");
  Int_t ncomb=0;
  
  const TGraphErrors *gtmp = g2;
  for(Int_t ii=0; ii<gtmp->GetN(); ii++){
    Double_t xx=-999, yy=-999;
    gtmp->GetPoint(ii,xx, yy);
    if(xx> -0.5 && xx< 0.5){
      Double_t ex = 0;//gtmp->GetErrorX(ii);
      Double_t ey = gtmp->GetErrorY(ii);
      
      hcomb->SetPoint(ncomb, xx, yy);
      hcomb->SetPointError(ncomb, ex, ey);
      printf("point %d: %f %f\n", ncomb, xx, yy);
      ncomb++;
    }
  }

  gtmp = g1;
  for(Int_t ii=0; ii<gtmp->GetN(); ii++){
    Double_t xx=-999, yy=-999;
    gtmp->GetPoint(ii,xx, yy);
    if(xx>0.5 && xx<2.7){
      Double_t ex = 0;
      if(gtmp->GetErrorX(ii)>0){
        ex = gtmp->GetErrorX(ii);
      }
      Double_t ey = gtmp->GetErrorY(ii);
      
      hcomb->SetPoint(ncomb, xx, yy);
      hcomb->SetPointError(ncomb, ex, ey);
      printf("point %d: %f %f\n", ncomb, xx, yy);
      ncomb++;
    }
  }
  
  gtmp = opt.Contains("pp")?g1:g2;
  for(Int_t ii=0; ii<gtmp->GetN(); ii++){
    Double_t xx=-999, yy=-999;
    gtmp->GetPoint(ii,xx, yy);
    if(xx> 2.7 && xx< 3.5){
      Double_t ex = 0;
      Double_t ey = gtmp->GetErrorY(ii);
      
      hcomb->SetPoint(ncomb, xx, yy);
      hcomb->SetPointError(ncomb, ex, ey);
      printf("point %d: %f %f\n", ncomb, xx, yy);
      ncomb++;
    }
  }

  hcomb->Draw("ap");

  Double_t fPar[]={0.208205, 1.859973, 7.089334, 0.141273 , 6.460552,  0.005147, -1.629302,  0.890983 };
  const Int_t npar = 8;
  
  TF1 *ff=new TF1("ff",AliTRDdEdxUtils::MeandEdxTRLogx, -0.5, 5, npar);
  ff->SetParameters(fPar);
  ff->Draw("same");
    
  XGLUtils::ChisquareFit(hcomb->GetN(), hcomb->GetX(), hcomb->GetY(), 0x0, hcomb->GetEY(), AliTRDdEdxUtils::MeandEdxTRLogx, npar, fPar);
  
  TCanvas *cc=new TCanvas(Form("c%s",opt.Data()));
  hcomb->SetMarkerStyle(20);
  hcomb->SetMarkerSize(1);
  hcomb->Draw("ap");
    
  ff->SetParameters(fPar);
  ff->SetLineColor(kRed);
  ff->Draw("same");
    
  cout<<endl;
  printf("opt: %s\n", opt.Data());
  gSystem->Exec("date +%m%d%Y%H%M%S");
  for(Int_t ii=0; ii<npar; ii++){printf("par[%d]=%15e;\n", ii, ff->GetParameter(ii));} cout<<endl;cout<<endl;

  
  TFile *fout=new TFile(Form("fittrd%s.root",opt.Data()),"recreate");
  hcomb->Write();
  ff->Write();
  cc->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char *argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=3){
    printf("argc!=3\n"); return 1;
  }

  const TString fin(argv[1]);
  TFile::Open(fin);

  TString opt(argv[2]);
  if(fin.Contains("Q0"))
    opt+="Q0";
  else if(fin.Contains("Q1"))
    opt+="Q1";

  if(opt.Contains("tpc"))
    fittpc();
  else
    fittrd(opt);

  return 0;

}
